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A quantum Monte Carlo method of determining Jastrow-Slater wave functions for which the 
energy is stationary with respect to variations in the single-particle orbitals is presented. A potential 
is determined by a least-squares fitting of fluctuations in the energy with a linear combination of 
one-body operators. This potential is used in a self-consistent scheme for the orbitals whose solution 
ensures that the energy of the correlated wave function is stationary with respect to variations in 
the orbitals. The method is feasible for atoms, molecules, and solids and is demonstrated for the 
carbon and neon atoms. 

PACS numbers: 71.15.-m, 31.25.-v, 02.70.Lq 
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Over the past decade, quantum Monte Carlo (QMC) 
methods [Q-|| have been used to calculate the structural 
and electronic properties of a variety of atoms, molecules, 
clusters and solids. For systems with large numbers of 
electrons, QMC methods at present provide the most ac- 
curate benchmark calculations of structural energies. In 
both variational (VMC) and diffusion (DMC) calcula- 
tions, a key step is the construction of a trial correlated 
many-electron wave function. In many such calculations, 
the trial wave function '^ is chosen to be of the Jastrow- 
Slater form ||-||,|1, i.e., * == JD, where D is a de- 
terminant of single-particle orbitals and J is the (posi- 
tive) Jastrow correlation factor. Although considerable 
progress has been made (principally using the variance 
minimization approach) in the numerical construction of 
optimal Jastrow factors plR-P] , relatively little attention 
has been given to the physical understanding and numer- 
ical optimization of the antisymmetric part of the wave 
function. In few-electron systems, variance minimiza- 
tion has been applied to the optimization of the determi- 
nant [|7| but, in larger systems, local density functional 
(LDA) or Hartree-Fock (HF) orbitals have generally been 
used as the only practical and numerically accurate way 
of constructing the determinantal part of the wave func- 
tion [||-|j|. It has not been clear why LDA (or HF) 
orbitals, which have little formal justification in this con- 
text, do so well or how one might in practice do better. 
A physical understanding of this issue and a practical 
method of approach to the calculation of such orbitals is 
bound to be particularly important to a wide and suc- 
cessful application of QMC methods. 

In this Letter, a new iterative method is demon- 
strated which successively updates the determinant D 
in Jastrow-Slater wave functions so that, at convergence, 
the energy of the full correlated wave function is sta- 
tionary with respect to variations in the single-particle 
orbitals in the Slater determinant. The method is cast in 
the framework of a self-consistent field problem for the 
determinant and can make use of the standard numerical 
codes (either LDA or HF), combined with VMC sampling 



methods of many-body wave functions. 

We will first derive Euler-Lagrange equations satis- 
fied when the energy is stationary with respect to the 
single-particle orbitals in the determinant. Then, we will 
present the numerical approach for an iterative scheme 
which solves those equations exactly [[l^l , demonstrate it 
in some numerical examples, and show that it can be con- 
veniently combined with existing variance minimization 
methods for the optimization of the Jastrow factor. 
Euler-Lagrange Equations: We assume that S is held 
fixed and that only the single-particle orbitals (j)i of the 
determinant D are varied. A general infinitesimal varia- 
tion of the orbital 4>i can be written as 4>i -^ 4>i+^j Vji't'j 
where rjji are infinitesimal coefficients and the sum is over 



a set of orthonormal orbitals which excludes all 



al- 



ready in D. The corresponding variation in D is given by 
Z) — > [1 -|- ^j rjjiC;Ci\D where cj and Ci are the fermion 
creation and destruction operators for the orbitals j and 
i, respectively. The energy is stationary with respect to 
all variations of the orbitals (f>i if and only if the Euler- 
Lagrange equations. 






(1) 



are satisfied for all i and j. ?i is the many-body Hamil- 
tonian. Explicit evaluation of the derivatives shows that 
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where E = ($|H|^) and I^P) is assumed normalized. 

By considering arbitrary linear combinations of the 
variations r]ji, we can see that solving the set of Euler- 
Lagrange equations (0) is equivalent to requiring that 



AEo = {^\{H-E) 
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for all possible one-body operators O. As discussed be- 
low, it is sometimes convenient or sufficient to consider a 



restricted class of variations of the orbitals in D so that 
Eq. (||) is satisfied only for a restricted class of one-body 
operators. 

Iterative Solution of the Euler- Lagrange Equations: We 
wish to find D such that AEo^ = for a set of n 
one-body operators {Ok}- We sample Nc configurations 
{R{i)} with local energies {E{i) = 7^*(i?(i))/*(i?(i))} 
from the square of the wave function \1/ = JD. We per- 
form a least-squares fit of the local energies with the sum 
£^0 + J2k=i VkOk{i), where Eq and 14 are fitting param- 
eters and Ok{i) = OkD{R{i))/D{R{i)). In the limit of 
Nc — > oo, this is equivalent to minimizing the integral 



m\n-E,-Y^Vk^ 
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with respect to the fitting parameters. This is in turn 
equivalent to solving the set of linear equations: 



^ 14 {AOkAOi) = (ASAOz), for ; - 1, 



(5) 
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where (•) denotes the average over the N,. configurations, 
AOfe(i) = Ok{i) - (Ok) and AE{i) = E{i) - (E). For 
Nc — > oo, {AEAOi) — > AEo^ and the fitting coefficients 
Vk are all zero if and only if all AEq^ — 0. 

Suppose we determined a set of V^ for a wave func- 
tion "^^^^ — JD^^^ according to the above procedure. 
How do we use the coefficients V^. to obtain the orbitals 
(pi for the determinant D^'^^ of the next iteration? 

Let's first suppose we have a non-interacting system 
with Hamiltonian TicS and eigenfunctions i^i but that 
we start from the eigenstates (p\ of an incorrect Hamil- 
tonian H^^^ = Hcff — X]fe=i^fe^fe- If '^6 follow the 
above fitting procedure (Eq. H) and determine the co- 
efficients 14, it is easy to see that 14 = A^ for all k 
and Tieff rnay be found (if not known in advance) as 
n^e = W(i) + Sn where SH = XiLi VfcCfc. The cor- 
rect single-particle orbitals 4>j are then eigenstates of the 
Hamiltonian H^^^ + SH. 

Motivated by this argument for non-interacting sys- 
tems, we construct the determinant D^^' of the corre- 



lated wave function ^(^■' from orbitals 



a(i) 



which are 



eigenstates of a non-interacting Hamiltonian Ti^^^. We 
compute V^: from Eq. ra and, as in the non-interacting 

case, use SH^^^ — J2k=i ^fe ^k as an increment to the 
physical external potential Vext in Ti.^^' to determine a 
set of orbitals ipl ■ ^ new increment STi.^-^^ is sim- 
ilarly derived from J'D^^^ and the external potential 
Vext + STi'^^^ + 5TiS^^ is used for the next iteration to ob- 
tain the orbitals cpi . Convergence is reached when 61i}-^'> 
is negligible. For accelerated convergence, the orbitals 
at the mih iteration are determined by perform- 



lim) 



ing a standard self-consistent LDA calculation with the 



external potential 14a;t + X^I^T '^^^'■*- The scheme con- 
verges within two or three iterations for the applications 
studied here. It should be noted that the LDA potential 
here is used purely for computational convenience and 
that the final orbitals (pf , the Hartrce potential and the 
final energy-fluctuation potential (EFP) part of the non- 
interacting Hamiltonian, H^'^p = ^™ ^ JH^'' + Vjf"^, 
are independent of the LDA. 

One may wonder why this approach, based on an argu- 
ment for non-interacting systems, would give rapid con- 
vergence in the interacting case. To see this, we con- 
sider the combined action of the many-body Hamilto- 
nian H and the Jastrow factor J7 on the determinant D 
by defining HjD{R) = [HJ{R)D{R)]/ J{R). The eigen- 
values of TL and Tij are identical and the eigenfunctions 
^ of H satisfy $ = J f where / is the corresponding 
eigenfunction of Hj. With a suitable choice of J , the 
two-body terms in TLj can be made weak ]1^. For the 
uniform electron gas, Bohm and Pines determined the 
long-wavelength part of the two-body Jastrow factor to 
remove the long-range fiuctuations of the two-body in- 
teraction in Ttj pTj . Similarly, the short-range cusp con- 
dition |1^ removes the e^ jvij divergence in Jij. If Jij 
were truly a non-interacting Hamiltonian Tioff, its exact 
eigenfunctions would be Slater determinants D. Thus, 
the motivation for choosing a trial Jastrow-Slater wave 
function JD is directly related to the approximation that 
two- and higher-body terms in 'Hj can be neglected Q . 
This, in turn, motivates our iterative approach to the 
solution of the Euler-Lagrange equations (||). 

In order to specify the full numerical implementation of 
the method, we need to choose an appropriate set of one- 
body operators Ok- We consider three possible choices: 

(1) Ok is a local potential /fc(r) so that Okii) = 
Sjli/fe(i"j) foi' tl^e configuration R{i) = (ri, . . . ,rAr). 
The evaluation of such an operator and of the aver- 
ages required in Eq. ^ and the increment of 14a;t by 
W{v) = X]fe=i ^fe/fe(r) in the self-consistent LDA cal- 
culation is then straight-forward. This case corresponds 
to the variational freedom of multiplying all orbitals (\)i 
in Z3 by a common function 1 -|- r]kfkij) and is equivalent 
to minimizing the energy with respect to the one-body 
term in the Jastrow factor p^ . 

(2) Ok is an angular-momentum-dependent potential 
fkifyPi with Vi the projection operator for angular mo- 
mentum I- The evaluation of the coefficients in Eq. @ can 
be made using the standard methods for the integration 
of non-local pseudopotentials in VMC H] and increment- 
ing the external potential in the LDA code by an angular- 
momentum-dependent potential is also straight-forward. 
This case corresponds to the variational freedom of mul- 
tiplying different angular momentum orbitals (\)i in D by 
different factors. 

(3) Arbitrary variations of the orbitals 0; may be al- 
lowed by using the one-body operators Oji — CjCi, where 



i labels an occupied orbital of the determinant D and j 
an unoccupied orbital. Oji acting on D simply replaces 
the orbital 0, in D with the orbital (f)j. The averages 
in Eq. [Scan be efhciently calculated using the relations 
in Ref. ]13| for replacing a row in a Slater determinant. 
The Hamiltonian in the self-consistent LDA calculation 
is then incremented by dTi — ^ ■ ■ Vji\(j)j){(f)i\ + c.c. 

Thus, while the approach can use, with trivial modi- 
fications, all the computational techniques to determine 
LDA orbitals, the final orbitals 4>i minimize the en- 
ergy for the many-body wave function JD^''^' with no 
restriction on the form of J . 
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FIG. 1. The incremental EFP potentials, 5H{r), for the 
ground state of the carbon and neon pseudo-atoms. The po- 
tentials corresponding to the Jastrow factors J\ and J2 are 
shown for the first and last iteration. 87i{r) at the first iter- 
ation is calculated using LDA orbitals in the determinant. 

In this Letter, we present results for the carbon and 
neon pseudo-atoms (pseudopotentials are used to elim- 
inate the Is core electrons 13]). The application of the 
approach to extended systems will be discussed in detail 
elsewhere ||lj]. In optimizing the orbitals, we have in- 
vestigated method (1) with a local EFP potential, and 
method (2) with both s and p EFP potentials. Because 
there is only one type of s orbital (spin up and down) and 
only one radial function for the p orbitals, methods (2) 
and (3) are equivalent and all variations of the orbitals 
consistent with the ground state symmetry are allowed 
in method (2). For simplicity, we will focus on results 
obtained using the local potentials (1). Interestingly, 
the orbitals, charge density, and energies differ only very 
slightly when the full variational freedom of the orbitals 
is allowed using separate s and p non-local potentials. 

The initial atomic orbitals in the Jastrow-Slater wave 
function are determined from a LDA calculation. Be- 
cause of self-interaction in the LDA, 6H{r) computed at 
the first iteration must behave like — e^/r at large dis- 
tances. Since, at large radii, the sampling of S7i has 
large statistical noise due to the very low electron den- 



sity, we constrain (57Y*^^^(r) to behave like — e^/r at large r 
while allowing full variational freedom at smaller r. This 
is achieved by writing SH'-^'^ (r) = Vo{r) + J^Hi Vk fk{r), 
where /fc(r) = cos[(/c — l)7rr/rc] exp[— (r/r^)''] and Vo(r) 
goes like — e^/r for r > r^ and smoothly becomes con- 
stant for r < Tc- The parameters Vk are determined by 
least-square fitting, as in Eq. H. After the first iteration. 
Voir) is not included in fitting 6'H^^'' . 

We performed the calculations for two different types 
of Jastrow factor, Jx and J2. The Jastrow factor 
Ji only contains electron-nucleus and electron-electron 
terms (see Appendix A of Ref. [|5[) and the value of 
its single free parameter was determined by minimiz- 
ing the energy. J2 includes electron-electron, electron- 
electron-nucleus and electron-nucleus terms (modified 
from Ref. M to deal with a pseudo-atom) and variance 
minimization was used to optimize its 25 parameters. 

In Fig. |l|, we show the first and last incremental EFP 
potentials for carbon and neon, obtained using the two 
Jastrow factors Ji and J2- The cut-off radius re is equal 
to 3 a.u. for neon and 5 a.u. for carbon and the number of 
basis functions, n/, is always less than 50. In each case, 
the final iteration is almost indistinguishable from zero, 
except for statistical sampling noise near the origin. The 
potential depends on the choice of the correlated compo- 
nent of the wave function: the superiority of J2 over J7i 
is reflected in the much smaller initial potential STL^^^r). 

TABLE L Total energies in VMC (Evmc) and DMC 
(Edmc) and root mean square fiuctuation [a) of the local 
energy in VMC. E^ and E^ are the percentages of correla- 
tion energies recovered in VMC and DMC. Hartree units are 
used. (For carbon, the HF energy is Shf = —5.3530 and the 
ground state energy |19] is Eq — —5.4561 Hartree. For neon, 
Shp = -34.6930 and Eq = -35.0106 Hartree pi). 
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-34.9554(2) 
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82.6(1) 
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JTi-Defp 


-34.9674(3) 
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86.4(1) 
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0.882 




J2£'lDA 


-34.9912(2) 


-35.0041(2) 


93.9(1) 


98.0(1) 


0.630 




J2DFFP 


-34.9913(2) 


-35.0040(2) 


93.9(1) 


97.9(1) 


0.624 



From Table I, we see that for each atom and for each 
type of Jastrow factor, the energy is lowered in going 
from LDA to EFP orbitals. Since J2 has greater flexibil- 
ity than J7i, this lowering of energy is negligible for ^2- 
For carbon (an open shell system) a multideterminant 
wave function is required to accurately represent the cor- 
relations and the use of a more flexible Jastrow factor J2 
gains little over J7i . (The very small difference in energy 
between ^2£'efp and Ji-Defp rnay be due either to the 



different parametric forms of J\ and J2 or to intrinsic 
differences between variance minimization and energy 
minimization.) In DMC, the energy gain in using EFP 
instead of LDA orbitals is negligible for both systems. 
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FIG. 2. Valence charge density 47rr^/9(r) for carbon and 
neon. The upper panel for each atom shows the "exact" den- 
sity p(|. The lower panel shows the differences from the HF 
density of the "exact", LDA, GGA (PW91) and the densi- 
ties for the wave functions ^i-Defp (VMCl) and ^2-Defp 
(VMC2). 

In Fig. g, the densities are shown for carbon and neon. 
In both systems and for both Jastrow factors, the VMC 
density obtained with the EFP approach is substan- 
tially closer to the best estimate of the true density than 
HF, and much closer than either LDA or the Perdew- 
Wang '91 (PW91) generalized gradient approximation 
(GGA) [|r^. In carbon, neither J\ nor J2 can capture 
the intrinsic multi-configuration correlation and the ac- 
curacy in the densities cannot rival that obtained in neon 
with i72Defp- 

In conclusion, we have demonstrated, for the first 
time, a numerically stable, rapidly convergent method 
which combines Monte Carlo sampling with existing self- 
consistent field techniques to minimize the energy with 
respect to the orbitals in a correlated Jastrow-Slater wave 



function. The approach may be combined with variance 
minimization methods for the optimization of the Jas- 
trow factor. The resulting variational many-body wave 
functions have electron densities very close to the most 
accurate densities available for atoms, using variance 
minimization and DMC methods. Preliminary tests show 
that the approach, in modified form, is also applicable to 
multi-determinant wave functions. We thank C. Umri- 
gar for useful discussions and E. Shirley for the use of his 
Hartree-Fock code. This work was supported by Enter- 
prise Ireland, Contract SC/98/748. 
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